function [optfunc,optfunc_frz] = polfunc_function(param)

        %=========================================================================%
        % Part 1: parameters
        %=========================================================================%
        rng(1,'twister');
        s = rng;
 	    sigma = param(1);      %EIS      
	    beta  = .97;           %discount factor 
        aT    = 81;            %terminal age
        aS    = 50;            %start age
        r     = .029;          %interest rate 
        pi    = .028;          %inflation        
        V_term = 0;            %terminal value
        z_age_min = 56;        %earliest freeze age
        z_age_max = 64;        %latest freeze age
        N_obs  = 5000;         %Number of simulated individuals
        ages = aS:1:aT-1;      %ages
        t = (1:(aT-aS))';      %time periods
        scaler = 1;            %scaler on consumption utility 
        solver_param = [sigma;beta;aT;aS;r;V_term;scaler;pi]; %collect some params
      
        %-------------------------------------------------------------------------%
        %Calculate annuity value and earnings per period from HRS
        %-------------------------------------------------------------------------%
        %1. HRS db wealth by age
        pdata = csvread('compwealth_frz_w10_22.csv',1,0);
        %col2 is age
        %col3 is avg db wealth in 2010 dollars
        %col4 is avg earnings in 2010 dollars
        dbw = pdata(aS-48+1:aT-48,2:3);
        y_raw = pdata(aS-48+1:end,4);
        %Smooth the earnings profile using cubic polynomic fit
        y_fit_param =  polyfit(ages',y_raw,3);
        y = polyval(y_fit_param,ages');
        %percent earnings losses from freezes (periods 0-4)
        y_cut = [.038;.021;.015;.049;.016];
        %2. SS annuity by retirement age (hh-level)
        b_ss = csvread('ssw10.csv',1,0);
        b_ss = b_ss(aS-48+1:end,4);
        %3. SSA mortality rates by age (48+)
        pdeath = csvread('SSA_mortality_2010.csv',1,0);
        pdeath = pdeath(aS+1:end,:);
        avpdeath = mean(pdeath(:,2:3),2);
        %4. annuity factors
        afactor = zeros(length(avpdeath),1); 
        %loop over age
        for a = 1:length(afactor)
            disc = zeros(length(afactor)-a+1,1);
            disc(1) = 1;
            for i = 2:length(disc)
                disc(i) = (1-avpdeath(a+i-1))*(1/((1+r)*(1+pi)))^(i-1); 
            end

            afactor(a) = sum(disc);
        end     
        afactor = [pdeath(1:aT-aS,1) afactor(1:aT-aS)];
        %DB annuity by retirement age
        b_db = dbw(:,2)./afactor(:,2);
        %total annuity benefit
        b = b_db + b_ss;
        
        %-------------------------------------------------------------------------%
        % Normalized mortality rate
        %-------------------------------------------------------------------------%
        mortp = avpdeath(1:aT-aS);
        mortp_80 = mortp./avpdeath(aT-aS+1);
        
        %-------------------------------------------------------------------------%
        % Asset grids (A = non-pension; W = DC pension)
        %-------------------------------------------------------------------------%
        %A
        A_gridN = 20;
        split = 5;
        A_mins = [0;50000;150000;500000;1000000;2000000];
        As=zeros(A_gridN/split,split);
        for i = 1:5
            step = (A_mins(i+1)-A_mins(i))/(A_gridN/split);
            As(:,i) = A_mins(i):step:(A_mins(i+1)-step);
        end
        A = reshape(As,[A_gridN,1]);
        
        %W
        W_gridN = 20;
        split = 5;
        W_mins = [0;25000;75000;250000;500000;1500000];
        Ws=zeros(W_gridN/split,split);
        for i = 1:5
            step = (W_mins(i+1)-W_mins(i))/(W_gridN/split);
            Ws(:,i) = W_mins(i):step:(W_mins(i+1)-step);
        end
        W = reshape(Ws,[W_gridN,1]);    
        
        %-------------------------------------------------------------------------%
        % Average tax rates from nber taxsim (age-dependent grids)
        %-------------------------------------------------------------------------%       
        worktax = csvread('workinc_taxsim_a_20.csv',0,0);
        %col1 = age, col2 = earn_before_dc, col3 =fica, col4=earn_after_dc
        %col5=fiitax, col6=dc_contrib
        itau_w = worktax(:,5)./worktax(:,4); %income tax
        ptau_w = 0.5*worktax(:,3)./worktax(:,2); %worker share of FICA
      
        rettax = csvread('retinc_taxsim_a_20.csv',0,0);
        %col1 = fiitax, col2 = age, col3=pension distribution, col4 = SS
        itau_r = rettax(:,1)./(rettax(:,3)+rettax(:,4));
             
        %Collect parameters for model solution
        taxf = cell(3,1);
        taxf{1,1}=itau_w;
        taxf{2,1}=ptau_w;
        taxf{3,1}=itau_r;
        
        %freeze-affected average tax rates 
        taxfz = cell(9,3);
        zindex = 1;
        for aa = 56:64         
            i_z = find(ages==aa);
            
            %While working
            fname = strcat('workinc_taxsim_a_20_z',num2str(i_z),'.csv');
            worktax_z = csvread(fname,0,0);
            itau_w_z = worktax_z(:,5)./worktax_z(:,4); %income tax
            ptau_w_z = 0.5*worktax_z(:,3)./worktax_z(:,2); %worker share of FICA
            taxfz{zindex,1} = itau_w_z;
            taxfz{zindex,2} = ptau_w_z;
     
            %Retirement
            fname = strcat('retinc_taxsim_a_20_z',num2str(i_z),'.csv');
            rettax_z = csvread(fname,0,0);
            itau_r_z = rettax_z(:,1)./(rettax_z(:,3)+rettax_z(:,4));
            taxfz{zindex,3} = itau_r_z;
            zindex = zindex+1;
        end
        
        %-------------------------------------------------------------------------%
        %Disutility of work process 
        %-------------------------------------------------------------------------%
    	mu      = 0;  %mean normal underlying f
        rho     = 1/(1+exp(param(2)));  %persistence of random process 
    	sigma_v = exp(param(3)); %sd of normal underlying f
    	phi     = param(4); %age slope
        gamma   = param(5); %constant
    	f_N     = 6; %grid size

        %=========================================================================%
        % Part 2: solve model for value and policy function 
        % without any freeze activity
        %=========================================================================%

        %Cell array to store value and policy functions for each age
        optfunc = cell(length(t),8);

        %-------------------------------------------------------------------------%
        %Solve age T-1 policy and value functions outside the loop
        %-------------------------------------------------------------------------%	
        %Set last period=true
        last=1;
        a=0;
        [V_w,A_polw,W_polw]=vwork_5(A,W,last,a,y,optfunc,solver_param,...
                    f_N,mu,rho,sigma_v,phi,gamma,taxf,mortp_80);
        [V_r,A_polr,W_polr]=vret_4(A,W,last,a,b_db,b_ss,optfunc,solver_param,taxf,mortp_80);
	
        
        % Outer value function and labor supply choice
        V = zeros(size(V_w)); 
        L = zeros(size(V_w)); 
       
        for i = 1:length(A)
            for j = 1:length(W)
                for k = 1:f_N
                    V(i,j,k) = max(V_w(i,j,k),V_r{aT-aS-a}(i,j));
                    L(i,j,k) = (V_w(i,j,k)>V_r{aT-aS-a}(i,j));
                end
            end
        end

        %Store value and policy functions
        optfunc(1,:) = {V_w;A_polw;W_polw;V_r;A_polr;W_polr;V;L};
        
        %-------------------------------------------------------------------------%
        %Solve policy and value functions for remaning periods
        %-------------------------------------------------------------------------%      
        %Set last period = false
        last=0;
        
        %Loop through ages
        for a = 1:length(t)-1 
    

               [V_w,A_polw,W_polw]=vwork_5(A,W,last,a,y,optfunc,solver_param,...
                    f_N,mu,rho,sigma_v,phi,gamma,taxf,mortp_80);
               [V_r,A_polr,W_polr]=vret_4(A,W,last,a,b_db,b_ss,optfunc,solver_param,taxf,mortp_80);


		% Outer value function and labor supply choice
		V = zeros(size(V_w)); 
		L = zeros(size(V_w)); 

		for i = 1:length(A)
		    for j = 1:length(W)
			for k = 1:f_N
			    V(i,j,k) = max(V_w(i,j,k),V_r{aT-aS-a}(i,j));
			    L(i,j,k) = (V_w(i,j,k)>V_r{aT-aS-a}(i,j));
			end
		    end
		end

		%Store value and policy functions
        optfunc(a+1,:) = {V_w;A_polw;W_polw;V_r;A_polr;W_polr;V;L};

        end        

        %=========================================================================%
        % Part 3: solve model for value and policy function 
        % separately for each potential freeze age 
        %=========================================================================%	
        %Cell array to store value and policy functions for each freeze age
        % dimension is age x freeze ages x number of functions (8)
        optfunc_frz = cell(length(t),z_age_max-z_age_min+1,8);
        lastfuncW0 = [];
        lastfuncR0 = [];
        
        parfor z = 1:z_age_max-z_age_min+1 %%loop over possible freeze ages
            
            %Freeze the path of DB pension annuities (keep SS wealth fixed)
            i_z = find(ages==z_age_min-1+z);
            b_db_z = b_db;
            b_db_z(i_z:end) = b_db(i_z);
            %freeze fixes nominal annuity; compute real annuity at each
            %prospective retirement age
            for i = i_z:length(b_db_z) 
                b_db_z(i) = b_db_z(i)*(1/(1+pi))^(i-i_z);
            end
            b_z = b_db_z + b_ss;
            %earn profile given freeze at age z
            y_loss = ones(length(y),1);
            y_loss(i_z:i_z+length(y_cut)-1) = 1-y_cut;
            y_z = y.*y_loss;

            %---------------------------------------------------------------------%
            %Solve age T-1 policy and value functions outside the age loop
            %---------------------------------------------------------------------%
            %Set last period=true
            last=1;
            a=0;

            [V_w,A_polw,W_polw]=vworkz_5(A,W,last,a,y_z,lastfuncW0,solver_param,...
                        f_N,mu,rho,sigma_v,phi,gamma,taxfz,mortp_80,z);
            [V_r,A_polr,W_polr]=vretz_4(A,W,last,a,b_db_z,b_ss,lastfuncR0,solver_param,taxfz,mortp_80,z);


            % Outer value function and labor supply choice
            V = zeros(size(V_w)); 
            L = zeros(size(V_w)); 

            for i = 1:length(A)
                for j = 1:length(W)
                    for k = 1:f_N
                        V(i,j,k) = max(V_w(i,j,k),V_r{aT-aS-a}(i,j));
                        L(i,j,k) = (V_w(i,j,k)>V_r{aT-aS-a}(i,j));
                    end
                end
            end

            %Store value and policy functions
            funcs = {V_w,A_polw,W_polw,V_r,A_polr,W_polr,V,L};   

            %-------------------------------------------------------------------------%
            %Solve policy and value functions for remaning periods
            %-------------------------------------------------------------------------%
            %Set last period = false
            last=0;

            %Loop through ages
            for a = 1:length(t)-1 

                lastfuncW = cell2mat(funcs(a,7));
                lastfuncR = funcs(a,4);

                [V_w,A_polw,W_polw]=vworkz_5(A,W,last,a,y_z,lastfuncW,solver_param,...
                        f_N,mu,rho,sigma_v,phi,gamma,taxfz,mortp_80,z);
                [V_r,A_polr,W_polr]=vretz_4(A,W,last,a,b_db_z,b_ss,lastfuncR,solver_param,taxfz,mortp_80,z);


                % Outer value function and labor supply choice
                V = zeros(size(V_w)); 
                L = zeros(size(V_w)); 

                for i = 1:length(A)
                    for j = 1:length(W)
                    for k = 1:f_N
                        V(i,j,k) = max(V_w(i,j,k),V_r{aT-aS-a}(i,j));
                        L(i,j,k) = (V_w(i,j,k)>V_r{aT-aS-a}(i,j));
                    end
                    end
                end

                %Store value and policy functions
                funcs0 = {V_w,A_polw,W_polw,V_r,A_polr,W_polr,V,L};
                funcs = [funcs;funcs0];

            end

            optfunc_frz(:,z,:) = funcs;

        end
               
end      

